function TSMC_Line_Correlate
%#ok<*NASGU>
%#ok<*ASGLU>
%#ok<*UNRCH>

GoldenInputImage='C:\Users\Les Frackowiak\Desktop\Tests\TSMC\Samples\46-G.jpg';
SampleInputImage='C:\Users\Les Frackowiak\Desktop\Tests\TSMC\Samples\46-S.jpg';
UseBackgroundEqualization=1;
BackgroundBlurSize=15;
BackgroundBlurSigma=10;
UseGrayscaleEqualization=1;
AbsoluteGrayscaleValue=150;
TrimAmount=30;
Method=1;
% Method = 0 : Mean Value
% Method = 1 : STD
% Method = 2 : Thresholded Sum
GrayThresholdLevels=[0 128];
% GrayThresholdLevels=[Minimum Maximum]
FilterType=0;
% FilterType = 0 : None
% FilterType = 1 : Median
% FilterType = 2 : Savitzky-Golay FIR
FilterSize=11;
% Must be odd integer
FilterOrder=3;
% Must be less than FilterSize-1 for Savitzky-Golay FIR Filter
NumberOfOutputs=10;


%%%% IMAGE LOADING %%%%
GoldenImage=imread(GoldenInputImage);
SampleImage=imread(SampleInputImage);
if size(size(GoldenImage),2)==3
    GoldenGray=rgb2gray(GoldenImage);
else
    GoldenGray=GoldenImage;
end
if size(size(SampleImage),2)==3
    SampleGray=rgb2gray(SampleImage);
else
    SampleGray=SampleImage;
end
[GoldenHeight,GoldenWidth]=size(GoldenGray);
[SampleHeight,SampleWidth]=size(SampleGray);
MinWidth=min(GoldenWidth,SampleWidth);
MinHeight=min(GoldenHeight,SampleHeight);
GoldenGray=GoldenGray(1:MinHeight,1:MinWidth);
SampleGray=SampleGray(1:MinHeight,1:MinWidth);
GoldenGray=GoldenGray(TrimAmount:(MinHeight-TrimAmount),TrimAmount:(MinWidth-TrimAmount));
SampleGray=SampleGray(TrimAmount:(MinHeight-TrimAmount),TrimAmount:(MinWidth-TrimAmount));

%%%% IMAGE PREP %%%%
if UseBackgroundEqualization
    Filter=fspecial('gaussian',BackgroundBlurSize,BackgroundBlurSigma);
    GI_Filtered=imfilter(GoldenGray,Filter,'symmetric');
    SI_Filtered=imfilter(SampleGray,Filter,'symmetric');
    GI_Filtered_Min=min(min(GI_Filtered));
    SI_Filtered_Min=min(min(SI_Filtered));
    GoldenGray=GoldenGray-(GI_Filtered-GI_Filtered_Min);
    SampleGray=SampleGray-(SI_Filtered-SI_Filtered_Min);
end
if UseGrayscaleEqualization
    GI_AverageValue=mean2(GoldenGray);
    SI_AverageValue=mean2(SampleGray);
    GoldenGray=GoldenGray*(AbsoluteGrayscaleValue/GI_AverageValue);
    SampleGray=SampleGray*(AbsoluteGrayscaleValue/SI_AverageValue);
end

%%%% 1D CONVERSION %%%%
if Method==0
    GoldenHorizontal=mean(GoldenGray,1);
    SampleHorizontal=mean(SampleGray,1);
    GoldenVertical=mean(GoldenGray,2);
    SampleVertical=mean(SampleGray,2);
elseif Method==1
    GoldenHorizontal=std(double(GoldenGray),0,1);
    SampleHorizontal=std(double(SampleGray),0,1);
    GoldenVertical=std(double(GoldenGray),0,2);
    SampleVertical=std(double(SampleGray),0,2);
elseif Method==2
    GoldenThresholded=im2bw(GoldenGray,GrayThresholdLevels(1)/255)-im2bw(GoldenGray,GrayThresholdLevels(2)/255);
    SampleThresholded=im2bw(SampleGray,GrayThresholdLevels(1)/255)-im2bw(SampleGray,GrayThresholdLevels(2)/255);
    GoldenHorizontal=sum(GoldenThresholded,1);
    SampleHorizontal=sum(SampleThresholded,1);
    GoldenVertical=sum(GoldenThresholded,2);
    SampleVertical=sum(SampleThresholded,2);
end
if FilterType==1
    GoldenHorizontal=medfilt1(GoldenHorizontal,FilterSize);
    SampleHorizontal=medfilt1(SampleHorizontal,FilterSize);
    GoldenVertical=medfilt1(GoldenVertical,FilterSize);
    SampleVertical=medfilt1(SampleVertical,FilterSize);
elseif FilterType==2
    GoldenHorizontal=sgolayfilt(GoldenHorizontal,FilterOrder,FilterSize);
    SampleHorizontal=sgolayfilt(SampleHorizontal,FilterOrder,FilterSize);
    GoldenVertical=sgolayfilt(GoldenVertical,FilterOrder,FilterSize);
    SampleVertical=sgolayfilt(SampleVertical,FilterOrder,FilterSize);
end
Horizontal=xcov(GoldenHorizontal,SampleHorizontal,'coef');
Vertical=xcov(GoldenVertical,SampleVertical,'coef');
SortedHorizontal=sort(Horizontal,'descend');
SortedVertical=sort(Vertical,'descend');
Horizontals=SortedHorizontal(1:NumberOfOutputs);
Verticals=SortedVertical(1:NumberOfOutputs);
H_Indices=zeros(NumberOfOutputs,1);
V_Indices=zeros(NumberOfOutputs,1);
for k=1:NumberOfOutputs
    H_Indices(k)=find(Horizontal==Horizontals(k),1,'first');
    V_Indices(k)=find(Vertical==Verticals(k),1,'first');
end
H_Shift(:,1)=[1:NumberOfOutputs];
V_Shift(:,1)=[1:NumberOfOutputs];
H_Shift(:,2)=H_Indices-length(GoldenHorizontal);
V_Shift(:,2)=V_Indices-length(GoldenVertical);
H_Shift(:,3)=Horizontals;
V_Shift(:,3)=Verticals;
H_Shift
V_Shift
end